////////////////////////////////////////////////////////////////////////////////////////////////////////////
//                                                                                                        //
//  ----------------------------- Ebrahim Foulaadvand, 12 May 2012 -----------------------------------    //
//                                                                                                        //
//     The programme "PopulationGrowth" evaluates the temporal evolution of a population number which     //
//     is subjected to stochastic linear birth and death processes.                                              //
//                                                                                                        //
//                                                                                                        //
////////////////////////////////////////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <conio>
#include <math.h>
#include <numeric>
#include <values>
#include <vector>
#include <list>
#include <algorithm>
#include <numeric>
#include<iomanip.h>

using namespace std;


double  tau=0.001,a=1,b=0.04; // tau=time step, a=birth rate, b=death rate

int N0=5,T=10000,t; // N0=initial number of individuals,

main()
{

ofstream  fileN ("N-t a=1 b=0.04 N0=5.plt"); //output file for computed time dependecne of population number.
ofstream  fileNanal ("Nanal-t a=1 b=0.04 N0=5.plt"); //output file for analytic time dependecne of population number.


vector<double> N(T+1,0),Nanal(T+1,0);

// Array N stores the computed number of individuals at time step t.
// Array Nanal stores the analytic number of individuals at time step t.


N[0]=N0;  // ---------------- initial conditions --------------------


//---------------------- Euler-Cromer Method    ---------------------

for(int t=0;t<T;t++){

N[t+1]=N[t]*(1+a*tau-b*tau*N[t]); //cout<<"N[t]= "<<N[t]<<endl;getch();

fileN<<t*tau<<" "<<N[t]<<endl;

fileNanal<<t*tau<<" "<<N0*a*exp(a*t*tau)/(a+b*N0*(exp(a*t*tau)-1))<<endl;

}

cout<<"end."<<endl;getch();

return 0;

}